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Preface 


This  research  grew  out  of  a  desire  to  create  a  highly  reflective  multilayer  film  with 
very  low  resistance  for  vertical  cavity  surface  emitting  semiconductor  lasers.  The 
traditional  cavity  mirrors  for  these  devices  are  formed  by  alternating  layers  of  gallium 
arsenide  and  aluminum  gallium  arsenide.  The  abrupt  changes  in  material  which  make  a 
good  mirror  also  make  a  poor  conductor,  so  electrical  pumping  of  these  lasers  is  difficult. 
One  solution  is  to  design  a  mirror  with  no  abrupt  changes  in  material  (lowering  the 
resistance)  while  maintaining  the  high  reflectance  required  for  the  laser  cavity  mirrors.  In 
the  course  of  the  investigation,  the  original  goal  of  designing,  fabricating,  and  testing  a 
vertical  cavity  laser  device  was  found  to  be  infeasible  in  the  time  available.  My  research 
therefore  took  a  more  theoretical  path,  and  I  decided  to  address  the  more  general  problem 
of  gradient  index  thin  film  design. 

I  would  like  to  express  my  thanks  to  my  advisor,  Maj  Jeffiey  Grantham,  for  his  aid 
and  support  in  my  ever  evolving  dissertation.  Dr  Peter  Haaland  was  also  of  great  help 
during  the  development  of  the  SWIFT  algorithm.  I  would  also  like  to  thank  Maj  Gregory 
Warhola  for  his  willingness  to  adopt  me  as  my  dissertation  topic  took  a  more 
mathematical  turn.  Finally,  I  wish  to  thank  my  wife  Joan  for  her  patience,  understanding, 
and  support  throughout  the  years  it  took  to  complete  my  studies. 
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Abstract 


Gradient  index  thin  films  provide  greater  flexibility  for  the  design  of  optical  coat¬ 
ings  than  the  more  conventional  “layer”  films.  In  addition,  gradient  index  films  have 
higher  damage  thresholds  and  better  adhesion  properties.  In  this  dissertation  I  present  an 
enhancement  to  the  existing  inverse  Fourier  transform  gradient  index  design  method,  and 
develop  a  new  optimal  design  method  for  gradient  index  films  using  a  generalized  Fourier 
series  approach.  The  inverse  Fourier  transform  method  is  modified  to  include  use  of  the 
phase  of  the  index  profile  as  a  variable  in  rugate  filter  design.  Use  of  an  optimal  phase 
fiinction  in  Fourier-based  filter  designs  reduces  the  product  of  index  contrast  and  thickness 
for  desired  reflectance  spectra.  The  shape  of  the  reflectance  spectrum  is  recovered  with 
greater  fidelity  by  suppression  of  Gibbs  oscillations  and  shifting  of  side-lobes  into  desired 
wavelength  regions.  A  new  method  of  gradient  index  thin  film  design  using  generalized 
Fourier  series  extends  the  domain  of  problems  for  which  gradient  index  solutions  can  be 
found.  The  method  is  analogous  to  existing  techniques  for  layer  based  coating  design,  but 
adds  the  flexibility  of  gradient  index  films.  A  subset  of  the  coefficients  of  a  generalized 
Fourier  series  representation  of  the  gradient  index  of  refraction  profile  are  used  as 
variables  in  a  nonlinear  constrained  optimization  formulation.  The  optimal  values  of  the 
design  coefficients  are  determined  using  a  sequential  quadratic  programming  algorithm. 
This  method  is  particularly  well  suited  for  the  design  of  coatings  for  laser  applications, 
where  only  a  few  widely  separated  wavelength  requirements  exist.  The  generalized 
Fourier  series  method  is  extended  to  determined  the  minimum  film  thickness  needed,  as 
well  as  the  index  of  refraction  profile  for  the  optimal  film. 


Xll 


1.  Introduction 


1.1.  Motivation 

The  Air  Force  is  interested  in  optical  thin  film  coatings  for  a  variety  of  uses, 
including  broad  band  anti-reflection  coatings  for  low-observable  applications,  narrow  pass 
filters  for  sensor  protection,  and  high  power  laser  mirrors.  Most  current  thin  film  coatings 
are  designed  and  constructed  using  a  number  of  slabs  of  materials  with  different  indices  of 
refraction  and  thicknesses.  The  most  popular  method  consists  of  alternating  layers  of  high 
and  low  index  material.  The  properties  of  the  film  are  determined  by  the  two  index  values 
used  and  the  thicknesses  of  each  layer  in  the  film.  Recently,  there  has  been  a  growing 
interest  in  inhomogeneous  or  gradient  index  thin  films,  which  are  characterized  by  a 
continuously  varying  index  of  refraction  throughout  the  film.  The  advantages  of  such 
gradient  index  films  over  the  alternating  stack  type  films  are  based  on  the  elimination  of 
the  abrupt  interfaces  between  materials.  In  a  slab  based  design,  very  large  electric  fields 
can  develop  at  the  material  interfaces.  This  can  lead  to  damage  of  the  film  in  high  power 
laser  applications  [47:272].  Also,  the  current  deposition  methods  demonstrate  columnar 
growth  patterns  in  each  material.  The  columnar  growth  patterns  create  voids  in  the 
structure,  which  can  increase  scattering.  In  addition,  the  different  material  properties  lead 
to  large  inherent  stress  between  layers  in  the  film  [47:271-272].  Elimination  of  the 
interfaces  by  using  gradient  index  designs  promotes  better  adhesion  to  the  substrate, 
reduces  internal  stress  in  the  film,  reduces  scattering,  and  increases  the  damage  threshold 
for  the  film  [43:2864]. 

Several  basic  approaches  to  thin  film  design  exist,  including  graphical,  analytical, 
and  digital  design  methods  [11:97].  Digital  techniques  are  of  particular  interest  because 
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they  can  be  used  to  design  films  with  more  complex  properties  than  is  possible  analytically 
or  graphically.  Digital  design  techniques  can  be  categorized  as  either  refinement  or 
synthesis  techniques.  Refinement  techniques  begin  with  an  initial  film  design  and 
iteratively  refine  it  to  achieve  a  better  design.  Synthesis  techniques  on  the  other  hand 
generate  a  film  design  based  on  the  desired  film  characteristics  only.  Several  methods 
exist  for  both  refinement  of  an  initial  design  [11]  or  synthesis  of  a  film  [12].  One  of  the 
synthesis  techniques  relies  on  the  inverse  Fourier  transform  to  relate  the  index  of 
refraction  of  the  film  to  the  transmittance  of  the  film.  This  technique,  originally  developed 
by  Sossi  and  Kard  [34,35,36],  will  be  expanded  upon  in  this  work.  In  addition,  new 
design  techniques,  based  on  a  generalized  Fourier  series  approach,  will  be  developed. 

The  objective  in  designing  a  thin  film  is  to  take  the  performance  requirements  for 
the  film  and  generate  a  design  subject  to  constraints  of  available  indices  of  refi-action  and 
acceptable  total  thickness.  The  performance  requirements  are  usually  stated  as 
reflectance,  transmittance,  or  absorption  versus  wavelength,  incident  polarization  or  angle. 
The  conventional  variables  in  the  design  are  the  number  of  layers  in  the  film,  and  the  index 
and/or  thickness  of  each  layer.  For  "simple"  filters,  such  as  notch  or  bandpass  filters,  anti¬ 
reflection  coatings,  and  reflectors,  designs  by  graphical  methods  or  by  knowledge  of  the 
properties  of  periodic  multilayers  are  possible  [12,  24].  Films  with  more  complex 
properties  require  digital  design  methods. 


1.2.  Problem  Statement 

The  purpose  of  this  research  is  to  create  a  design  method  to  map  between 
reflectivity  and  continuous  index  of  refraction  profiles  that  also  allows  for  additional 
design  constraints,  such  as  limitations  on  the  total  film  thickness  and  index  of  refi’action 
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range.  The  process  of  designing  an  optical  filter  requires  a  mapping  between  the  space  of 
all  possible  index  profiles  and  the  space  of  all  possible  reflectivity  characteristics.  A 
mapping  from  index  to  reflectivity  can  be  generated  from  Maxwell's  equations.  This 
mapping  allows  one  to  evaluate  the  performance  characteristics  of  a  film,  in  terms  of 
reflectivity  as  a  function  of  wavelength,  angle  of  incidence,  polarization,  etc.. 
Unfortunately,  this  mapping  is  not  invertable,  so  it  cannot  be  used  to  design  a  film  (except 
by  trial  and  error).  Most  of  the  existing  design  techniques  rely  on  the  trial  and  error 
approach.  The  corrections  to  the  index  are  based  on  numerical  principles  of  optimal  design 
theory.  One  exception  to  this  approach  is  the  inverse  Fourier  transform  technique.  This 
design  method  is  derived  from  Maxwell’s  equations  by  making  several  approximations. 
However,  this  method  maps  one  reflectivity  profile  to  many  possible  index  profiles.  In 
addition,  there  is  no  mechanism  to  impose  additional  constraints  on  the  design,  such  as 
available  index  range  or  total  thickness. 

This  research  is  presented  in  two  phases.  The  first  phase  is  to  investigate  design  of 
gradient  index  thin  films  by  the  inverse  Fourier  transform  method.  The  second  phase  is  to 
create  a  new  method  for  solving  the  inverse  problem  of  identifying  an  index  of  refi'action 
profile  for  a  given  reflectivity,  based  on  a  generalized  Fourier  series  approach.  Two 
variations  on  this  theme  are  explored:  a  Fourier  series  method  and  a  wavelet  method.  The 
wavelets  provide  a  different  framework  for  analyzing  the  structure  of  the  index  of 
refraction.  A  wavelet  decomposition  of  a  “signal”  is  analogous  to  a  Fourier  series 
decomposition,  but  the  information  is  packaged  in  terms  of  “scale”  and  “shift”,  rather  than 
the  more  familiar  frequency  and  phase.  Another  reason  to  use  wavelets  is  that  each 
wavelet  has  a  finite  extent,  as  opposed  to  the  sines  and  cosines  of  Fourier  series.  This 
finite  support  is  helpful  because  it  allows  one  to  focus  only  on  the  elements  of  the  design 
that  are  important.  As  an  example,  the  inverse  Fourier  transform  method  requires  the  user 
to  specify  the  reflectivity  over  a  broader  range  than  the  design  requirements  specify.  This 
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is  required  to  achieve  sufficient  detail  in  the  index  of  refraction  profile.  This  forces  the 
design  to  meet  more  stringent  conditions  than  are  really  needed.  The  finite  support  of  the 
wavelets  should  allow  one  to  overcome  such  problems.  The  theory  and  concepts  of  this 
wavelet  formalism  are  presented  in  more  detail  in  Appendix  C. 

1.3.  Organization 

Chapter  2  presents  background  information  on  existing  thin  film  design  techniques 
and  the  basic  theory  of  thin  films.  Chapter  3  discusses  the  inverse  Fourier  transform 
method  for  synthesis  of  thin  films,  as  well  as  a  modification  to  the  existing  theory  allowing 
optimal  synthesis  of  a  thin  film.  Chapter  4  discusses  optimal  design  of  gradient  index  films 
using  Fourier  series  and  wavelet  series.  The  final  chapter  summarizes  the  dissertation  and 
presents  suggestions  for  future  research  in  this  area.  The  appendices  include  an 
explanation  of  the  numerical  considerations  involved  in  implementing  the  inverse  Fourier 
transform  techniques  of  Chapter  3  (Appendix  A),  a  brief  outline  of  optimization  theory 
(Appendix  B),  the  wavelet  theory  needed  in  this  research  (Appendix  C),  and  the  programs 
used  in  this  work  (Appendix  D). 
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2.  Background 


This  chapter  presents  the  historical  background  of  thin  film  design  methods  and  the 
current  state  of  the  art  for  gradient  index  film  design  and  fabrication.  In  addition,  the  basic 
theory  of  optical  properties  of  thin  films  necessary  for  the  analysis  and  synthesis  of  thin 
films  is  presented. 


2.1.  Historical  Perspective 

To  appreciate  the  significance  of  the  new  approach  to  gradient  index  thin  film 
design  presented  in  this  dissertation,  a  brief  history  of  thin  film  design  techniques  is 
needed.  Thin  films  were  first  discovered  in  the  late  1600’s  by  Robert  Hooke  and  Sir  Isaac 
Newton  in  the  phenomenon  known  as  “Newton’s  rings”  [19:299,24:2].  The  first  anti¬ 
reflection  coatings  were  made  by  Fraunhaufer  in  1817  [24:2].  The  modem  era  of  thin  film 
manufacturing  began  in  the  1930’s,  with  the  invention  of  reliable  vapor  deposition 
techniques  [24:4]. 

Modem  design  techniques  can  be  broken  into  three  broad  categories;  analytical 
methods,  graphical  methods,  and  digital  designs.  Analytical  techniques  use  a  few  simple 
elements  as  building  blocks  to  design  more  complex  filters.  The  building  blocks  used,  such 
as  quarter  wave  stacks,  are  selected  because  the  relationship  between  the  few  variables  of 
the  basic  element  (i.e.,  index,  thickness,  number  of  layers)  and  the  reflectivity  (i.e.,  total 
reflectivity,  pass  band  width,  stop  band  location)  are  well  known.  A  specific  design  is  then 
created  by  combining  a  series  of  these  building  blocks  until  the  desired  profile  is  obtained 
[24:164-172].  This  technique  is  effective  for  fairly  simple  designs,  such  as  a  bandpass 
filter,  but  does  not  work  as  well  for  more  complex  designs.  Many  of  the  existing  thin  film 
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design  techniques  rely  on  the  intuition  of  the  designer,  based  on  his  experience  with  basic 
elements  of  the  design.  One  example  of  this  is  the  design  of  bandpass  filters  from  quarter 
wave  stacks  [24:164-172],  Another  example  is  the  design  of  rugate  filters,  where  the 
basic  building  blocks  are  sinusoids  [22],  It  has  been  noted  that  one  of  the  standard  rugate 
filter  profiles  for  a  bandpass  filter  bears  a  remarkable  resemblance  to  a  Morlet  wavelet 
[37]. 

Graphical  techniques  are  closely  related  to  the  analytical  methods  described  above. 
They  consist  of  various  methods  to  simplify  the  calculations  necessary  in  the  analytical 
method  by  graphing  the  relationships  between  the  variables.  Examples  of  these  techniques 
include  the  Smith  chart,  reflection  circles,  and  admittance  loci  [24:54-66]. 

The  techniques  of  primary  interest  here  are  the  digital  design  methods.  Digital 
design  methods  can  be  further  subdivided  into  two  categories;  refinement  and  synthesis 
techniques.  Refinement  techniques  are  characterized  by  the  fact  that  they  require  an  initial 
starting  design.  Synthesis  techniques,  on  the  other  hand,  create  their  design  without  an 
initial  guess  [11:  2876].  There  are  over  a  dozen  published  refinement  methods,  and  about 
half  a  dozen  synthesis  methods.  There  are  undoubtedly  many  more  unpublished  design 
methods  that  are  proprietary  to  the  thin  film  manufacturing  companies.  Most  of  the 
methods  make  use  of  a  merit  function,  which  is  a  measure  of  how  close  the  existing  design 
is  to  the  desired  result.  Some  of  the  methods  are  capable  of  locating  a  global  minimum, 
which  is  the  true  optimal  solution  to  the  problem.  Other  methods  only  solve  for  local 
minima  based  on  the  initial  guess.  Brief  descriptions  of  a  few  of  the  published  design 
methods  are  given  below. 

Refinement  Methods:  [11:2878-2882] 

1.  Adaptive  Random  Search.  This  procedure  takes  a  starting  design  and  applies  a 
random  change  to  the  parameters  of  the  design.  The  results  at  each  step  are  compared  to 
the  desired  result  through  a  merit  function.  Changes  in  the  design  that  reduce  the  merit 
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function  are  kept,  while  others  are  discarded.  The  magnitude  of  the  change  in 
parameters  is  also  changed  at  each  step,  making  this  procedure  perform  both  a  local  and 
global  search.  This  means  that  if  the  “best”  design  is  very  different  from  the  starting 
design,  this  method  may  find  it.  Local  search  methods  only  make  small  changes  to  the 
initial  design.  This  method  stops  its  iteration  when  the  merit  function  passes  a  preset 
threshold  or  the  method  stops  converging 

2.  Damped  Least  Squares:  This  method  uses  the  derivatives  of  the  quantities  of 
interest  in  the  merit  function  with  respect  to  the  construction  parameters  to  determine  the 
changes  to  be  made  to  the  design.  Since  the  design  problem  is  highly  non-linear,  the  size 
of  the  change  in  any  one  step  is  limited  by  a  damped  least  squares  algorithm.  This  method 
performs  a  local  search  from  the  starting  design. 

3.  Golden  Section  Method:  The  golden  section  method  takes  and  initial  design  and 
varies  each  of  the  construction  parameters  in  sequence.  The  optimal  value  for  each 
parameter  is  found,  (within  the  preset  limits  for  each)  before  continuing  on  to  change  the 
next  parameter.  This  algorithm  converges  fairly  rapidly,  and  can  find  solutions  far  from 
the  initial  design. 

4.  Hookes  and  Jeeves  Pattern  Search:  This  technique  is  composed  of  two  steps; 
an  exploration  and  a  pattern  search.  The  exploration  is  to  change  each  parameter  up  and 
down  by  a  small  amount,  and  determine  which  changes  improve  the  merit  function.  The 
pattern  search  extrapolates  in  the  direction  of  improvement  from  the  exploration  step. 
The  step  sizes  are  changed,  and  then  the  process  is  repeated.  The  process  stops  when  the 
change  in  the  merit  function  is  smaller  than  a  preset  limit. 

5.  Simulated  Annealing:  The  physical  process  of  annealing  is  heating  a  sample  and 
allowing  the  molecules  to  find  a  minimum  energy  state  during  cooling.  This  idea  is  used 
mathematically  in  this  design  approach.  This  method  can  deal  with  arbitrarily  large 
numbers  of  design  parameters,  and  can  find  globally  optimal  solutions. 
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Synthesis  Methods:  [23:3790-3791] 

1 .  Comprehensive  Search  Method.  This  method  is  limited  to  films  with  only  a  few 
layers  (<6).  The  global  optimal  solution  is  found  by  testing  every  possible  combination  of 
index  and  thickness  allowed.  Clearly,  only  a  limited  number  of  indices  and  thickness  can 
be  used.  The  advantage  of  this  method  is  that  it  finds  the  truly  optimal  solution  given  the 
constraints  on  index  and  thickness.  The  disadvantage  is  it  takes  so  much  computer  time  to 
check  every  possible  combination  that  only  very  simple  systems  can  be  designed  this  way. 

2.  Gradual-Evolution  Method:  This  method  is  a  modification  of  the 
comprehensive  search  described  above.  Instead  of  designing  the  whole  film  in  one  step,  a 
smaller  design,  say  three  layers,  is  done  by  comprehensive  search.  The  result  is  then  added 
as  part  of  the  substrate,  and  an  additional  few  layers  are  designed  on  this  “new”  substrate, 
again  by  comprehensive  search.  Since  only  a  few  layers  are  changed  at  a  time,  this  method 
is  much  faster  than  the  full  comprehensive  search  method.  However,  there  is  no  longer  a 
guarantee  that  the  global  optimal  solution  will  be  found. 

3.  Minus  Filter  Method:  A  minus  filter  is  a  design  that  transmits  all  incident 
radiation  except  that  in  a  narrow  spectral  band.  Thelen  has  shovm  that  all  dielectric  minus 
filters  can  be  made  using  quarter  wave  stacks  with  several  indices  of  refraction  [39:365- 
369].  This  design  technique  breaks  the  desired  reflectivity  profile  into  several  minus 
filters,  and  then  places  the  individual  designs  in  series. 

4.  Flip-Flop  Method:  This  unique  approach  uses  many  thin  layers  and  only  two 
indices  of  refraction.  The  algorithm  is  to  pass  through  the  film  in  sequences,  flipping  each 
layer  between  high  and  low  index  and  calculating  the  resultant  merit  function.  After 
several  passes  through  the  film,  this  method  converges  on  a  solution.  The  advantage  to 
this  technique  is  it  uses  only  two  materials,  which  is  preferred  by  manufacturers,  and  it 
requires  no  initial  design.  The  method  sometimes  results  in  alternating  very  thin  layers, 
but  this  can  usually  be  corrected  manually  later. 
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2.2.  Gradient  Index  Film  Design  and  Fabrication 


The  current  gradient  index  film  design  techniques  focus  on  the  inverse  Fourier 
transform  method  and  the  rugate  film  design  method.  The  inverse  Fourier  transform 
method  can  be  classified  as  a  synthesis  technique,  since  no  initial  design  is  required.  The 
rugate  design  method  is  an  analytical  design  technique. 

1 .  Inverse  Fourier  Transform  Method.  The  use  of  Fourier  transforms  in  the 
synthesis  of  thin  films  was  first  proposed  by  Delano  in  1967  [9].  Sossi  is  generally 
credited  with  the  development  of  the  inverse  Fourier  transform  method  for  thin  film  design 
[34,35,36].  Sossi's  papers  introduce  a  Fourier  transform  relationship  between  the 
logarithmic  derivative  of  index  of  refraction  and  the  reflectance  of  the  film.  This  design 
technique  was  further  developed  by  Dobrowolski  and  Lowe  in  1978  [10].  Many  authors 
have  applied  this  method  to  various  design  problems,  and  offered  modifications  to  the 
theory  to  better  suit  their  needs.  These  include  design  of  wideband  anti-reflectance 
coatings  [44]  and  high  reflectance  filters  [16,43],  Other  work  has  focused  on  modifying 
the  theory  of  the  method  to  improve  the  quality  of  the  results  or  ease  of  computation 
[4,6,14,42].  This  method,  along  with  a  modification  of  my  design,  will  be  described  in 
detail  in  Section  3.2  below. 

2.  Rugate  Design  Method.  A  rugate  filter  is  composed  of  a  number  of  basic 
elements  combined  to  achieve  the  desired  optical  properties  of  the  film.  The  fundamental 
rugate  design  element  is  a  sine  wave  refractive  index  profile  on  a  substrate  [22],  The  basic 
rugate  film  has  a  notch  reflectance  profile,  with  the  location,  width,  and  height  of  the 
notch  controlled  by  five  parameters.  More  complicated  rugate  films  are  created  by  adding 
several  basic  rugates  in  parallel  or  in  series,  and  including  apodization  functions  or 
matching  layers  to  suppress  sidebands  in  the  reflectance. 
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The  fabrication  of  films  designed  using  these  techniques  is  done  in  one  of  two 
ways.  The  first  way  is  to  use  the  concept  of  the  Herpin  equivalent  index  to  convert  a 
gradient  index  design  into  a  two  index  material  system,  and  make  an  equivalent  film  out  of 
discrete  layers  [3,15,24:191-200],  This  method  allows  the  use  of  gradient  index  design 
methods  to  design  conventional  alternating  layer  based  films,  but  does  not  offer  the 
advantages  an  actual  gradient  index  film  would.  The  other  fabrication  possibility  is  to 
deposit  the  gradient  index  profile  as  designed.  This  is  currently  on  the  forefront  of 
existing  technology,  and  is  an  active  area  of  research.  There  are  several  material  systems 
that  can  be  used  to  produce  gradient  index  profiles.  Table  2-1  lists  some  of  the  materials 
available,  along  with  the  index  range  and  literature  reference  for  additional  information. 


Table  2-1 :  Gradient  Index  Film  Materials 


Material  System 

Index  Range 

Reference 

Gei-xSx 

3.5 -4.0 

Ali-xGaxAs 

2.9 -3.6 

1-x^^x 

.  2.2 -2.5 

[22:97] 

SiOxNv 

1.5 -2.0 

[22:97,  29:179] 

Si02  -  Ta205 

1.65-2.0 

MgF2  -  Ti02 

1.38-2.3 

CeFs  -  ZnS 

1.6 -2.2 

1.38-2.5 

[20:197-204] 
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2.3 .  Optical  Properties  of  Thin  Films 

This  section  presents  the  theory  and  methods  for  determining  the  reflectance 
properties  of  a  film.  The  first  section  presents  the  derivation  of  the  reflectance  fi'om 
Maxwell’s  equations,  and  the  second  section  outlines  the  implementation  of  this 
formulation  for  a  single  thin  film.  The  third  section  extends  this  implementation  to 
gradient  index  films. 

2.3.1.  Derivation  of  Reflectance 

This  section  presents  the  formal  derivation  for  determining  the  reflectance  of  a  film 
starting  from  Maxwell’s  equations.  Figure  2.1  shows  the  coordinates  used  in  the 
derivation.  The  surface  of  the  film  is  the  x-y  plane,  and  the  physical  thickness  of  the  film  is 
along  z.  The  origin  is  at  the  film-substrate  interface,  and  the  thickness  of  the  film  is 
denoted  by  L. 


Figure  2. 1 :  Geometry  of  Thin  Film 
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Assuming  a  time  dependence  of  exp[/ffl/].  Maxwell’s  equations  in  SI  units  are 


curl  H  =  icosE 
curl  E  =  -ico/jH 
div  £  =  0 
div  H^O 


(2.1) 


where  E  is  the  electric  field,  H  is  the  magnetic  field,  s  is  the  permittivity,  and  //  is  the 
permeability. 

Two  polarizations  must  be  considered  to  determine  the  electric  field;  the 
transverse  electric  (TE)  case,  which  has  its  electric  field  aligned  perpendicular  to  the  plane 
of  incidence,  and  the  transverse  magnetic  (TM)  case,  which  has  its  electric  field  aligned 
parallel  to  the  plane  of  incidence.  By  the  symmetry  of  Maxwell’s  equations,  the  TM  case 
can  be  deduced  fi-om  the  TE  case  with  a  substitution  of  E  for  H  and  s  for  //.  The  plane  of 
incidence  is  defined  by  the  direction  of  propagation  of  the  incident  light  and  the  normal  to 
the  surface.  Define  the  x  direction  as  the  direction  in  which  the  electric  field  is  aligned,  so 
the  fields  can  be  written  as 


Ej,  (z)  exp[/(s;?  - 

0 

E  = 

0 

H  = 

Hy  (z)  exp[icot  -  //:5y] 

0 

(z)  exp[icvt  -  ikSy] 

Here  S  is  an  invariant,  S  =  nQ  sin(0o),  where  no  is  the  index  of  refraction  of  the  incident 
medium,  and  do  is  the  angle  of  incidence.  Also,  co  is  the  frequency  of  the  incident  light,  k 
is  the  wave  number  of  the  light,  k=27r/?i,  and  X  is  the  wavelength  in  vacuum.  Substituting 
these  expressions  for  E  and  H  into  Maxwell’s  equations  yields 
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dz 

dz 

H,{z)^ 


-  -ik 


^2\ 


MrJ 


=  -ikfi^ZoHJz) 


EAz) 


(2.3) 


where  Zq  is  the  impedance  of  free  space,  Zq  = 

^0.  Ma  are  the  permitivity  and  permeability  of  free  space, 

£r>  Mr  are  the  relative  permittivity  and  permeability  of  the  medium. 

Let  N^(z)  =  £^(z)jUr  -  S^,  and  consider  only  dielectric  materials  with  //r=l.  N(z) 
is  always  a  real  valued  function  since  5”  <  1  for  all  cases  considered  here  and  £>  >  1  for  all 
dielectric  materials.  For  normal  incidence,  S=0,  and  N  reduces  to  the  usual  definition  of 
index  of  refraction.  Also  define  two  new  complex  valued  variables  u  and  v  by 


u{z)  = 


v{z)  =  -MrZo-Y~ 


(2.4) 


The  normalization  factor  Et  is  the  complex  amplitude  of  the  transmitted  wave.  Using  the 
new  normalized  field  variables.  Maxwell’s  equations  reduce  to  the  following  pair  of 
coupled  first  order  ordinary  differential  equations 

u’{z)=ikv{z)  and  v’{z)-ikN^{z)u{z)  (2.5) 

with  boundary  conditions 


w(0)  =  1  v(0)  =  cos(^,) 


(2.6) 
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In  this  formulation,  the  amplitude  reflectance  and  transmittance  of  the  film,  r(k) 
and  t(k)  respectively,  are  found  from  the  Fresnel  equations:  [40:4274] 

/ ,  k)  cos(^o )  +  ^{L,  k)  cos(^, ) 

w(Z,  k)  cos((9o )  -  v(Z,  k)  cos(0, ) 

(2.7) 

, , .  ^  ln^,XL,k)cos{eo) 

k)  cos(^o  )  “  cos(^, ) 

Here  the  field  variable  dependence  on  the  wave  number  k  has  been  called  out  explicitly. 

2.3.2.  Multilayer  Thin  Films 

The  optical  properties  of  a  multilayer  thin  film  are  based  upon  the  Fresnel 
reflection  at  interfaces  and  the  interference  between  multiple  reflections.  The  reflectance 
of  such  a  film  can  be  determined  by  using  Maxwell's  equations  and  applying  the  boundary 
conditions  at  both  interfaces.  MacLeod  has  shown  the  fields  at  the  first  interface  are 
related  to  the  fields  at  the  second  interface  by  [24:35] 

cos(<^)  isin((^)/77i  ' ^2  8) 

HoiJ  l_i77isin(^)  cos{5) 

where  Eoi,  Hoi  represent  the  field  at  the  interface  between  the  incident  media  and  the  film, 
and  En,  Hn  represent  the  field  at  the  interface  between  film  and  the  substrate.  The  optical 
admittance  of  the  film  material  rj  is  given  by 
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TM 

cos{6) 

=  itiYo  cos{9)  TE 


(2.9) 


with  7o=1/377  Siemens.  The  phase  shift  <Jis  given  by 


S  =  Im^d  co^^)/l 


(2.10) 


where  rii  is  the  index  of  refraction  of  the  film,  d  is  the  thickness  of  the  film,  6  is  the  angle 
between  the  direction  of  propagation  and  the  normal  to  the  surface  of  the  film,  and  X  is  the 
wavelength  in  vacuum.  Note  that  the  phase  shift  5  can  be  complex,  and  the  angle  6^, 
found  using  Snell's  Law,  may  also  be  complex.  The  complex  amplitude  reflectivity,  p, 
and  the  real  intensity  reflectance,  R,  can  by  found  from  Equation  2.8  by  defining  an  input 
optical  admittance  for  the  assembly  as 


(2.11) 


Then  for  an  incident  medium  with  an  optical  admittance  of  t/q 


P  = 


rio-Y 

TJo+y 


rjo-Y'^ 

Ro+y) 


♦ 


(2.12) 


where  the  asterisk  indicates  complex  conjugation. 

Since  the  reflectance  depends  only  on  the  optical  admittance  of  the  film,  Y, 
rearrange  Equation  2.8  to  find  Y : 
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B 

cos(S) 

i  sin(<^)  / 

■  1  ■ 

C 

iTji  sin(<5) 

cos(S) 

.^2. 

Y  =  CIB 


(2.13) 


The  matrix  above  is  the  characteristic  matrix  for  the  film,  and  depends  only  on  the 
film  materials  and  the  angle  of  incidence  of  the  impinging  light.  This  matrix  technique  can 
be  directly  extended  to  an  assembly  of  many  thin  films.  So,  given  q  thin  film  layers,  each 
with  specified  index  and  thickness,  the  characteristic  matrix  of  the  assembly  is  given  by 


B 

C 


17] r  sin(^,.) 


i  sin(^,. )/  7]^ 
cos(Sr) 


(2.14) 


where  for  each  layer  r  the  phase  shift  Sr  is 

Sj.  =  2nnj.d^  cos(^^)/A 


(2.15) 


and  is  the  index  of  refraction  of  the  layer,  dr  is  the  thickness  of  the  layer,  and  the  angle 
Or  is  found  using  Snell's  Law.  These  equations  form  the  basis  for  the  design  of  thin  film 
dielectric  mirrors. 

In  addition  to  determining  the  reflectivity  of  a  multilayer  thin  film,  the 
characteristic  matrix  above  can  also  be  used  to  determine  the  phase  of  the  reflected  field. 
Writing  the  optical  admittance  of  the  film  as  Y=a  +  bi,  the  phase,  (p,  of  the  field  upon 
reflection  is  given  by  [24:37] 


(2.16) 
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The  simplest  example  of  a  mirror  using  these  multilayer  films  is  a  quarter  wave 
stack.  This  is  a  series  of  layers  each  with  an  optical  thickness  of  one  quarter  of  the 
wavelength  of  the  incident  light.  If  the  index  and  thickness  of  each  layer  are  chosen  so  the 
optical  path  length  nd  =  ^4,  and  the  light  is  at  normal  incidence,  the  phase  shift  5  for  each 
layer  is  Till.  The  characteristic  matrix  then  reduces  to 


B 

C 


n 

r=l  L' 


0  i/;?/ 

}Vr  0  , 


1 


(2.17) 


This  characteristic  matrix  is  very  easy  to  calculate.  An  example  of  such  a  quarter  wave 
stack  is  25  alternating  layers  of  GaAs  and  AlAs  on  a  GaAs  substrate,  which  have  indices 
of  refraction  of  3.6  and  3.2  respectively  [33:588].  If  the  design  wavelength  is  1.0  micron, 
this  mirror  is  93.6%  reflective  at  this  center  wavelength  and  has  a  reflectivity  versus 
wavelength  profile  as  shown  in  Figure  2.2. 


Figure  2.2:  Reflectivity  vs.  wavelength  for  25  layer  quarter  wave  stack  of 
GaAs/ AlAs  on  a  GaAs  substrate. 
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2.3.3.  Gradient  Index  Films 

The  analysis  above  assumes  the  mirror  is  formed  by  a  discrete  number  of  layers 
with  well-defined  thickness  and  index  of  refraction.  How  does  one  treat  a  structure  where 
the  index  varies  continuously?  One  approach,  derived  by  Bovard,  is  to  start  from  the 
beginning  with  Maxwell's  equations  and  develop  an  expression  for  the  characteristic 
matrix  of  the  film  based  on  the  known  index  of  refraction  as  a  function  of  position  in  the 
film  [4:1999-2001].  Unfortunately,  the  resultant  expression  contains  an  infinite  series  of 
integrals,  which  makes  evaluation  impracticable.  Another  alternative,  and  the  one 
implemented  here  (see  REFLECT.M  in  Appendix  D),  is  to  approximate  the  continuous 
film  as  a  discrete  film  with  many  small  layers.  As  a  general  rule,  discrete  layers  on  the 
order  of  5  nm  thick  give  a  good  approximation  for  the  reflectance  of  a  gradient  index  film 
for  visible  light  [6:5429].  Since  the  formalism  presented  above  includes  the  complex 
portion  of  the  fields  involved,  the  interference  effects  are  already  incorporated  in  the 
calculation.  This  is  of  critical  importance  to  the  validity  of  this  approach,  since  high 
reflectivity  depends  on  destructive  interference  of  the  transmitted  fields. 

The  reflectance  is  calculated  by  a  “C”  language  version  of  REFLECT.M.  The 
conversion  of  this  program  from  the  MATLAB™  language  to  “C”  increases  the  speed  of 
each  reflectance  calculation  by  at  least  a  factor  of  ten.  The  “C”  language  program, 
REFLECT.C,  is  also  included  in  Appendbt  D.  The  inputs  to  REFLECT.C  are  the 
sampled  values  of  the  index,  the  substrate  index,  the  vector  of  wavenumbers  k  at  which  to 
determine  the  reflectance,  and  the  array  of  sample  positions  in  the  film  x.  The  output  is  an 
array  of  reflectances,  one  for  each  of  the  input  wavenumbers  k.  The  units  for  x  and  k  must 
be  consistent,  i.e.  microns  and  inverse  microns  or  nanometers  and  inverse  nanometers. 
The  programs  REFLECT.M  and  REFLECT.C  are  used  throughout  the  dissertation  to 
determine  the  reflectance  characteristics  of  the  gradient  index  thin  film  designs. 
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Appendix  A  :  SWIFT  Numerical  Design  Considerations 


The  discrete  implementation  of  the  inverse  Fourier  transform  relation  in  Equation 
3.1  is  accomplished  using  the  fast  Fourier  transform  (FFT)  algorithm  in  the  MATLAB^“ 
programming  environment.  MATLAB^'^  is  a  “technical  computing  environment  for  high- 
performance  numeric  computation  and  visualization”  [27:/],  which  includes  a  high  level 
programming  capability.  The  MATLAB^^  programs  developed  for  this  work  are  included 
in  Appendix  D.  The  implementation  of  the  SWIFT  technique  is  done  using  two  programs: 
QUE.M  and  SWIFT.M.  The  program  QUE.M  uses  a  valid  expression  for  Q(k)  such  as 
one  of  the  forms  of  Equation  3.11  to  generate  a  sampled  version  of  Q(k)/k  for  notch 
reflectance  profiles.  SWIFT.M  calculates  the  phase  function  for  given  Q(k)/k  using  the 
SWIFT  technique  (Equations  3.13  and  3.18),  and  then  calculates  the  index  profile  using 
the  built  in  inverse  Fourier  transform  (Equation  3.1).  This  program  numerically  estimates 
the  integral  in  Equation  3.18  from  the  sampled  values  of  Q(k)  using  a  simple  trapezoid 
summation  algorithm.  A  third  program,  REFLECT.M,  is  used  to  estimate  the  reflectance 
for  a  gradient  index  film  by  approximating  the  film  by  thin,  homogeneous  layers  at  each 
sample  index  value  and  implementing  Equations  2.12  through  2.15. 

The  inputs  to  QUE.M  are  the  lower  and  upper  wavelength  limits  for  the  notch 
reflector  in  microns,  the  reflectivity  of  the  notch  (between  0  and  1),  the  number  of  samples 
to  use,  and  the  total  thickness  of  the  film  to  model  in  microns.  Note  that  the  total 
thickness  of  the  film  here  is  not  the  same  as  the  desired  film  thickness,  nor  the  “power 
spread”  thickness  of  the  SWIFT  theory.  The  relationship  among  these  three  design 
parameters  will  be  farther  discussed  below.  The  functional  dependence  of  Q(k)  on  the 
reflectivity  or  transmission  of  the  film  (such  as  in  Equation  3.11)  is  hard  coded  in  the 
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program.  The  output  of  the  QUE.M  program  is  a  vector  containing  the  samples  values  of 
Q(k)/k. 

The  inputs  to  the  SWIFT.M  program  are  the  vector  of  Q(k)/k  values  output  by 
QUE.M,  the  number  of  samples,  the  lower  and  upper  wavelength  values  in  microns  used 
to  generate  the  Q  ftinction,  the  total  film  thickness,  and  the  “power  spread”  thickness  for 
the  SWIFT  algorithm.  Again,  this  need  not  be  the  same  as  the  desired  film  thickness.  The 
output  of  the  program  is  a  matrix;  the  first  column  contains  the  sample  values  of  the 
index,  the  second  column  contains  the  optical  thickness  of  that  layer,  and  the  third  column 
is  the  position  of  the  sample  in  optical  thickness. 

The  inputs  to  the  REFLECT.M  program  are  the  index  matrix  of  index  values  and 
layer  thicknesses  from  SWIFT.M,  and  the  number  of  plot  points.  The  index  of  the 
substrate  and  exit  medium  are  fixed  in  the  code,  as  is  the  wavelength  range  to  calculate. 
The  output  is  a  matrix  with  the  reflectance  in  column  one  and  the  wavelength  in  microns  in 
column  two. 

There  are  a  few  numerical  computation  issues  associated  with  using  the  discrete 
implementation  of  the  Fourier  transform.  The  functions  QUE.M  and  SWIFT.M  above 
require  the  user  to  define  the  sampling  scheme  to  use  in  the  Fourier  transform  calculation. 
The  goal  is  to  accurately  specify  an  index  of  refi’action  profile  over  some  finite  thickness  of 
film.  There  are  several  related  parameters  that  must  be  specified  to  calculate  a  discrete 
Fourier  transform.  Note  that  the  Fourier  transform  variables  in  the  theory  are  the  double 
optical  thickness,  x,  of  the  index  of  refraction  and  the  wavenumber,  k,  of  the  Q  function. 
However,  the  reflectance  calculation  requires  the  index  profile  be  specified  in  terms  of 
optical  thickness.  For  consistency,  all  thicknesses  specified  as  inputs  to  and  outputs  from 
the  programs  are  in  optical  thickness,  and  the  conversion  to  double  optical  thickness  in 
done  inside  the  programs  as  necessary. 
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The  parameters  for  the  index  profile  are;  number  of  sample  values,  denoted  by  Ns, 
the  spacing  of  samples  in  optical  thickness,  denoted  by  A,  the  total  optical  thickness  of  the 
film  sampled,  drotai,  and  the  design  optical  thickness  of  the  film,  dfli„.  The  parameters  for 
the  Q  function  in  ^-space  are:  the  number  of  samples.  Ns,  the  spacing  of  samples,  6  = 
Nidrotai),  and  the  total  range  of  ^-space  sampled,  1/(A).  Clearly,  the  number  of  samples 
can  be  found  fi-om  the  simple  relation:  Ns  =  drotaJN  The  critical  frequency,  or  Nyquist 
frequency,  is  related  to  the  sample  spacing  of  the  index  by  /c=l/(2A).  The  various 
parameters  for  specifying  the  numerical  sampling  are  illustrated  in  Figure  A.  1  below.  The 
sampling  densities  actually  used  for  calculation  are  much  higher  than  depicted  in  this 
figure.  The  method  for  determining  the  sampling  density  needed  is  illustrated  in  the 
example  below.  All  the  examples  in  Section  3.3  use  the  same  approach. 


Figure  A.  1 :  Illustration  of  numerical  parameters  for  discrete  Fourier  transform. 
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Acceptable  values  for  the  design  parameters  are  found  by  considering  the  physical 
constraints  on  the  design  problem.  First,  the  number  of  samples,  Ns ,  should  be  a  power  of 
two  to  take  advantage  of  the  efficiency  of  the  fast  Fourier  transform  algorithm.  The 
sample  spacing  of  the  index  of  refraction  profile.  A,  should  be  on  the  order  of  one  sample 
every  five  nm,  to  ensure  the  reflectivity  calculation  is  accurate.  Also,  the  Q  function, 
which  is  the  starting  point  for  this  design,  must  be  adequately  sampled.  The  number  of 
non-zero  samples  of  the  Q  function  is  denoted  by  No.  Finally,  the  total  spatial  thickness  of 
the  film,  drotai,  must  be  sufficiently  large  to  minimize  spatial  aliasing  in  the  index.  Aliasing 
arises  from  the  finite  sampling  of  a  function.  Any  power  that  lies  outside  the  Nyquist  range 
is  folded  back  inside  the  range,  generating  in  this  case  an  inverse  Fourier  transform  that  is 
incorrect  [30:387].  Viewed  another  way,  this  means  that  the  sampling  of  the  Q  fiinction 
must  be  sufficiently  high  to  insure  the  index  profile  calculated  is  near  the  substrate  value  at 
the  edges  of  the  film. 

To  illustrate  the  process  to  determine  acceptable  parameter  values  for  a  problem, 
the  parameters  for  the  first  example  in  Section  3.3.1  will  be  found  below.  Recall  this 
example  is  to  design  a  narrow-band  reflectance  filter  with  a  reflectivity  of  90%  from  the 
lower  (initial)  design  wavelength,  X,^580,  to  the  upper  (final)  design  wavelength,  A,/=620 
nm,  and  0%  outside  this  band.  The  indices  of  the  substrate  and  the  incident  medium  are 
the  same  («,„&  =  Jiout  =  1.50).  The  desired  optical  thickness  of  the  film  is  df,im=30  pm.  The 
numerical  parameters  for  the  sampling  are  found  by  requiring  the  index  sample  spacing.  A, 
to  be  on  the  order  of  five  nm  in  optical  thickness,  the  number  of  non-zero  samples.  No,  of 
the  Q  function  to  be  about  25,  and  the  total  number  of  samples,  Ns ,  to  be  a  power  of  two. 
The  first  step  in  determining  appropriate  parameter  values  is  to  make  a  rough  estimate, 
based  on  the  values  above.  The  second  step  will  be  to  evaluate  the  numerical  values 
produced,  and  then  select  a  set  of  new  numerical  values  that  satisfy  all  the  constraints. 
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First,  since  the  sample  spacing  in  ^-space  is  the  inverse  of  the  total  thickness  (see 
Figure  A.l),  the  requirement  for  25  non-zero  samples  of  the  Q  function  over  notch  design 
wavelengths  can  be  used  to  determine  the  total  double  optical  thickness: 

N  25 

drotai^  X  =—l - j— =  225;mi  (A.l) 

Note  that  for  this  calculation  the  spatial  frequencies  are  defined  as  l/X,  rather  than 
2ti/X.  This  is  to  be  consistent  with  the  QUE.M  and  SWIFT. M  programs  in  Appendix  D. 
This  value  for  drotai  can  be  used  along  with  the  sample  spacing,  A,  of  the  index  profile  to 
estimate  the  number  of  samples  needed.  Ns.  The  sample  spacing.  A,  must  be  multiplied  by 
a  factor  of  2  to  convert  the  requirement  on  the  sample  spacing  in  the  index  from  optical 
thickness  (required  by  the  reflectance  calculation  program)  to  double  optical  thickness, 
which  is  the  variable  of  the  Fourier  transform  relation.  The  sample  spacing  value  to  use  in 
these  calculations  is  therefore  A=2(0.005)=0.010  pm. 


_  (A.2) 

At  this  point,  constrain  the  number  of  samples.  Ns  ,  to  be  a  power  of  two. 
Unfortunately,  the  Ns  found  above  is  not  a  power  of  two,  and  the  closest  powers  of  2  are 
2*'‘=16  384^  and  2^^32^786.  In  the  interest  of  minimizing  computation  time,  choose  the 
lower  number  of  samples,  and  then  determine  the  impact  on  the  other  parameters.  So,  for 
A^i=2*'‘=i6  384^  and  keeping  the  original  desired  sample  spacing  in  double  optical  path 
length,  A=0.010  pm,  the  new  total  film  thickness,  drotai,  would  be 
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=  NA 

=  (16,384)(0.010)  (A.3) 

=  163.84  nm 


Using  this  new  value  for  drotah,  the  number  of  non-zero  samples  of  the  Q  function  for  this 
parameter  choice  can  be  determined  by 


•^0  di^lgl 


j_ 

A, 


^  ■  ^0.58  0.6 


62 


18  non -zero  samples  of  Q(k) 


(A.4) 


At  this  point  the  designer  can  decide  to  accept  this  set  of  parameters,  or  decide  to  try  a 
different  combination  of  values.  For  example,  one  might  decide  that  18  samples  of  the  Q 
function  is  not  enough,  and  decide  to  adjust  the  total  thickness  again  to  increase  the 
number  of  samples.  This  would  also  affect  the  sample  spacing  of  the  index  profile.  To 
illustrate  this,  choose  to  keep  the  total  number  of  samples  Ai=2^'*=16,384,  but  now  choose 
to  increase  the  total  film  thickness  to  be  dTotai=  200  |am.  Using  this  new  value  for  djotai,  the 
number  of  non-zero  samples  of  the  Q  function  is 


■^0  dfg^gl 


j _ 


=  (200)1 


1 


1 


0.58  0.62/ 

23  non  -  zero  samples  of  Q{k), 


(A.5) 


and  the  sample  spacing  in  index,  A,  would  be 
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200 


16,384 

=  0.0122  //m 


(A.6) 


This  corresponds  to  an  optical  thickness  sample  spacing  of  0.0061  pm,  which  is  sufficient. 

This  collection  of  design  parameters  is  closer  to  the  original  objectives,  so  they  are 
selected  to  use  in  the  example.  The  final  parameters  required  as  input  to  the  programs  in 
Appendix  D  are  expressed  in  optical  path  length  (not  double  optical  path  length).  The 
conversion  from  optical  path  length  to  double  optical  path  length  is  done  in  the  programs 
as  needed.  The  reason  for  this  is  to  keep  the  input  consistent  with  that  required  by  the 
reflectance  calculation.  Therefore  the  final  parameter  values  for  total  film  thickness,  drotai, 
and  sample  spacing  in  index.  A,  are  the  values  found  above  divided  by  2.  For  the  example 
of  Section  3.3.1  the  parameters  are: 

2'^  =  16,384 

drotai  ^ '^00  fxm 

A  =  0.0061  nm 

Nq  =>  23  non  -  zero  samples  of  Q 

This  choice  for  the  total  optical  thickness  of  the  film  should  also  minimize  any  aliasing 
effects,  since  it  is  over  three  times  the  desired  film  optical  thickness.  The  same  process  is 
used  to  determine  the  parameters  used  in  all  of  the  examples  in  Chapter  3. 
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Appendix  B  :  Optimization  Theory 


Introduction 

The  process  of  design  requires  a  statement  of  the  problem  to  be  solved, 
identification  of  alternative  solutions,  and  some  measure  of  what  constitutes  a  “best” 
solution  to  the  problem.  For  simple  problems,  it  may  be  possible  to  exhaustively  list  all 
possible  solutions  and  “obviously”  determine  which  one  is  best.  For  more  complex 
problems,  however,  a  more  structured  approach  is  needed  to  ensure  the  “best”  solution  is 
found.  This  structured  approach  consists  of  building  a  mathematical  model  for  the  system 
in  question,  identifying  a  set  of  variables  to  adjust  in  the  design,  and  defining  a  “merit 
function”  to  numerically  identify  the  best,  or  “optimal”  solution.  This  section  will 
summarize  the  theory  of  optimal  design,  which  will  be  applied  to  the  problem  of  gradient 
index  thin  film  design  in  Chapter  4. 

Statement  of  Optimal  Design  Problem 

To  express  a  design  problem  in  mathematical  form,  one  first  must  define  the 
physical  system  to  be  modeled  and  identify  the  inputs  to  and  outputs  from  this  system. 
The  next  step  is  to  construct  a  mathematical  model  of  the  system  which  approximates  its 
performance  by  correctly  mapping  inputs  to  outputs.  The  model  consists  of  a  set  of 
design  variables  and  a  collection  of  functions  of  these  variables  that  define  the  objective 
and  the  constraints  on  the  problem.  The  objective  function  (sometimes  referred  to  as  the 
merit  function)  is  usually  constructed  such  that  its  minimum  value  represents  the  optimal 
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design.  The  formal  definition  of  the  constrained  optimization  problem  with  n  real-valued 
variables  is  therefore 


minimize  fix) 
jc  e9?” 

Subject  to: 

gi(x)  =  0  i  = 

gj(x)<0  i  = 


(B.l) 


where  jc  is  a  vector  of  the  design  variables, /(Oc^  is  the  objective  function,  the  gi(x)’s  are 
the  constraint  functions,  and  the  JCmin  and  jc^ax  are  vector  bounds  on  the  design  variables.  It 
should  be  noted  here  that  there  are  two  general  classes  of  constrained  optimization 
problems.  The  first  class  is  characterized  by  the  fact  that  all  the  functions  in  Equation  B.  1 
are  linear  functions  of  the  design  variables.  This  class  of  problems  can  be  readily  solved 
by  the  techniques  of  linear  programming.  The  second  class  of  problems,  and  the  one  with 
which  the  rest  of  this  paper  is  concerned,  have  non-linear  fijnctions  of  the  design  variables 
as  objective  functions  or  constraints.  This  format  of  the  optimal  design  problem  with  the 
constraints  all  expressed  as  less  than  or  equal  to  zero  is  called  the  Null  Negative  Form 
[28: 17].  The  solution  of  the  problem  expressed  by  Equation  B.  1  is  denoted  by  x*. 

A  few  additional  definitions/concepts  are  needed  before  describing  the  methods 
used  to  solve  the  optimal  design  problem.  First,  a  point  Xo  in  the  vector  space  9?”  is  a 
“feasible”  point  if  all  constraints  gfxo)  are  satisfied  by  equality  or  inequality,  as  required  by 
the  constraint.  Conversely,  a  point  in  91"  where  this  is  not  true  is  called  an  “infeasible” 
point.  The  second  concept  is  that  of  “active”  constraints.  A  constraint  is  said  to  be  active 
if  it  is  equal  to  zero  when  evaluated  at  the  solution  point  ( g,(x*)=0  ).  Clearly,  all  equality 
constraints  are  active  by  definition. 
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Now  that  the  mathematical  problem  has  been  defined,  how  is  a  solution 
determined?  Conceptually,  the  most  straightforward  method  would  be  to  select  a  feasible 
point  JCjtin  the  vector  space  SR”  as  a  start,  then  explore  the  neighborhood  by  adjusting  each 
element  of  the  vector  by  a  fixed  amount.  Then  evaluate  the  objective  function  and 
constraints  at  each  test  point,  and  identify  the  feasible  point  for  which  the  objective 
function  is  smallest.  This  test  point  would  become  the  new  starting  point,  Xk+i  ,  and  the 
procedure  would  be  repeated  until  a  minimum  value  for  the  objective  function  is  found. 
While  this  method  is  easily  visualized,  it  is  not  practical.  There  are  much  more  efficient 
techniques  for  determining  the  solution  to  the  constrained  optimal  design  problem.  The 
current  state  of  the  art  in  non-linear  optimization  is  a  method  known  as  Sequential 
Quadratic  Programming  (SQP)  [28:365-372].  The  basic  theory  for  this  method  is  outlined 
below. 

Kuhn-Tucker  Conditions  and  SQP 

The  general  approach  to  solving  non-linear  constrained  optimization  problems  is  to 
transform  the  problem  into  an  easier  sub-problem,  which  can  be  iterated  to  find  the  true 
solution.  The  foundation  of  this  approach  is  a  set  of  necessary  conditions  on  the  solution 
of  a  general  constrained  optimization  problem  known  as  the  Kuhn-Tucker  (KT)  equations: 
[26:2-22] 


m 

V/(**)  +  I/l  >,(:(•)=  0 
1=1 

i  =  (B.2) 

^  *>  0  i  =  m^  +  \,...,m 
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The  first  equation  is  the  gradient  of  the  Lagrangian,  and  describes  a  canceling  of  the 
gradients  between  the  objective  function,/,  and  the  constraints,  g,,  weighted  by  Lagrange 
multipliers,  Ai.  Only  active  constraints  are  included  in  this  expression,  so  the  Lagrange 
multiplier  for  any  inactive  constraint  is  set  to  zero.  The  solution  point  x*  is  assumed  to  be 
a  “regular”  point,  which  means  that  gradients  of  the  constraints  are  linearly  independent  at 
that  point.  Note  that  the  KT  equations  are  not  sufficient  conditions  for  a  solution,  so  each 
point  which  satisfies  Equation  B.2  must  be  checked  individually  to  determine  if  it  is  the 
optimal  solution.  An  additional  sufficiency  condition  can  be  formed  using  information  on 
the  second  derivative  of  the  Lagrangian.  For  any  regular  solution  point  of  the  KT 
equations,  x*,  if  the  Hessian  of  the  Lagrangian  at  jc*  is  positive-definite,  then  the  solution 
point  is  a  local  constrained  minimum  [28:184].  A  positive-definite  matrix  is  symmetric 
and  has  only  positive,  real  elements.  The  mathematical  expressions  for  the  Lagrangian 
(denoted  by  L)  and  the  Hessian  of  the  Lagrangian  (denoted  by  H),  are  given  below: 


Z(jc,  A  )  =  f(x)  +  i  ,.g,.(jc), 
1=1 


H{x,A)  = 


dx^dx^ 

dxjdx^ 


dx^dx^ 

d^L 

dxjdx^ 


dx^dx. 


d^L 

dx,dx„ 


d^L 

dx„dx„ 


(B.3) 


These  necessary  and  sufficient  conditions  for  the  solution  of  the  original 
constrained  optimization  problem  can  be  used  in  building  a  sequence  of  easier  sub¬ 
problems  to  solve. 
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The  easier  sub-problem  to  be  solved  is  to  determine  the  best  “direction”  in  the 
vector  space  of  variables  to  move  in  for  the  next  iteration,  as  well  as  the  length  of  step  to 
take  along  this  direction.  Again  using  the  null  negative  form  of  the  problem  in  Equation 
B.  1,  form  the  Lagrangian  as 


L{x,  1  )  =  f(x)  +  A  (x)  (B  .4) 

i=l 

The  simplified  problem,  called  the  Quadratic  Programming  (QP)  sub-problem,  is  to  solve 
the  quadratic  approximation  to  the  Lagrangian  above  subject  to  linearized  constraints 
[28:367],  The  quadratic  approximation  to  the  Lagrangian  is  found  by  expanding  in  a 
Taylor’s  series  and  keeping  terms  up  to  second  order  in  x.  For  any  fixed  point  in  the 
vector  space  of  variables,  x^ ,  the  QP  sub-problem  can  be  written  as 


minimize  +  'Vf(Xi^y  d 

subject  to: 

Vg,  (x^  fd  +  g,.  (x^ )  =  0,  /■  =  1,2, 

"^gXXkf  d  +  gXx^)<0,  i  =  m^  +l,...,/w 


Here  the  functions  /  and  g,  are  the  original  objective  function  and  constraints,  H  is  a 
positive-definite  approximation  to  the  Hessian  of  the  Lagrangian,  and  d  is  the  new  step 
direction  for  the  next  iteration.  Bold  typeface  indicates  a  vector  or  matrix.  The 
superscript  T  indicates  matrix  transpose.  The  solution  to  this  QP  sub-problem  is  used  to 
determine  the  next  iteration  point  =x^^+ad,  where  the  step  length  a  is  to  be 

determined.  The  step  length  a  is  needed  to  stabilize  the  algorithm,  because  it  is  possible 
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that  the  solution  to  the  QP  sub-problem  is  infeasible.  It  is  found  by  performing  a  line 
search  along  the  direction  d  and  minimizing  a  merit  function,  which  includes  penalties  for 
infeasible  solutions. 

In  summary,  the  SQP  algorithm  is; 

1.  Select  a  starting  point,  Xo,  and  an  initial  estimate  for  A,o.  (^0) 

2.  Solve  the  QP  sub-problem  above  to  determine  the  step  direction  d  and 
the  new  estimate  of  Xk+i-  (k=k+ 1 ) 

3.  SetXk^i  =  Xk+  ad. 

4.  Test  for  termination.  If  criteria  not  satisfied,  goto  step  2. 

The  SQP  algorithm  for  solution  of  non-linear  optimization  problems  is  available  in  the 
MATLAB^“  Optimization  Toolbox  [26,  27].  MATLAB'^^  is  an  interactive  software 
package  for  scientific  and  engineering  numeric  computation  [27:4].  The  next  section 
outlines  the  numerical  methods  used  by  MATLAB^“  to  implement  this  algorithm. 

MATLAB™  Implementation  of  SQP 

The  numerical  optimizations  performed  in  this  research  were  done  using  the 
MATLAB^*^  software  package  and  in  particular  the  MATLAB^*^  optimization  toolbox 
[26,27].  The  MATLAB^*^  software  uses  the  SQP  algorithm  described  above  in  its 
constrained  optimization  routines.  To  implement  the  theory  above,  the  program  must 
estimate  the  gradients  and  the  Hessian  at  each  point.  The  implementation  of  the  SQP  is 
done  in  three  main  stages  [26:2-26]: 

1 .  Update  the  estimate  for  the  Hessian  matrix  of  the  Lagrangian. 

2.  Solve  the  QP  sub-problem. 

3.  Perform  line  search  and  merit  function  minimization. 
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Equality  vs.  Inequality  Constraints 


The  MATLAB'T'^  algorithm  differentiates  between  equality  and  inequality 
constraints.  Equality  constraints  are  active  by  definition,  and  have  non-zero  Lagrange 
multipliers.  Inequality  constraints  are  included  in  null  negative  form,  and  may  or  may  not 
be  active.  The  program  uses  an  active  set  strategy  approach  to  handle  the  inequality 
constraints.  This  means  that  an  estimate  of  which  constraints  are  active  is  calculated  at 
each  iteration.  This  estimate  of  activity  figures  in  the  merit  function  used  in  the 
determination  of  step  size.  Also,  the  zero  tolerances  for  the  equality  and  inequality 
constraints  are  controlled  separately.  This  can  impact  the  results,  depending  on  how  the 
problem  is  originally  stated. 

Gradient  Estimation 

If  the  objective  functions  and  constraints  are  not  known  analytically,  the  partial 
derivatives  needed  to  estimate  the  gradients  must  be  calculated  numerically.  MATLAB'’’^ 
uses  an  adaptive  finite  difference  technique  to  find  the  values  of  the  gradients  as  needed. 
The  maximum  and  minimum  perturbations  for  this  calculation  can  be  controlled  by  the 
user. 

BEGS  Method  for  Hessian  Estimation 

The  key  to  the  SQP  algorithm  is  knowledge  of  the  Hessian  of  the  Lagrangian  for 
the  problem.  Since  this  matrix  is  computationally  expensive  to  calculate  directly,  it  is 
estimated  and  refined  through  iteration  using  the  formula  devised  by  Broyden,  Fletcher, 
Goldfarb,  and  Sharmo  (BFGS)  [26:2-26].  In  addition,  the  estimate  for  the  Hessian  is 
forced  to  remain  positive-definite,  so  the  solution  to  the  original  problem  is  a  minimum. 
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Termination  Criteria 


The  optimization  goal  is  achieved  when  three  termination  criteria  have  been  satis¬ 
fied.  There  are  three  numerical  tolerances  to  be  met,  one  on  the  objective  function,  one  on 
the  variables,  and  one  on  the  constraints.  The  objective  function  tolerance  specifies  the 
precision  required  on  the  objective  function  at  the  solution.  The  termination  criteria  for 
variables  specify  the  minimum  acceptable  precision  on  the  values  of  the  variables.  The 
constraint  termination  criterion  is  the  maximum  allowable  violation  of  the  constraints  at 
the  solution  point.  All  three  termination  criteria  must  be  satisfied  simultaneously  to 
achieve  an  optimal  solution.  The  default  termination  tolerances  are:  objective  ~  lO"'*  , 
variables  ~  10"* ,  and  constraints  —  10'^.  All  three  can  be  adjusted  if  necessary. 
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Appendix  C :  Wavelet  Theory 


Wavelets  are  a  relatively  new  mathematical  technique.  They  are  a  form  of 
generalized  Fourier  series,  using  wavelets  as  basis  elements  instead  of  sines  and  cosines. 
These  new  basis  elements  will  be  used  to  solve  Maxwell’s  equations  for  the  fields  inside 
the  film,  and  establish  a  new  form  for  the  mapping  between  index  of  refraction  and 
refiectivity.  The  sections  below  will  first  introduce  some  of  the  basic  concepts  of  wavelets, 
and  then  describe  the  technique  of  multi-resolution  analysis. 

Basic  Wavelet  Concepts 

The  term  wavelets  was  first  coined  in  1982  by  Morlet,  who  used  a  mathematical 
technique  similar  to  Fourier  series  but  using  “little  waves”  as  the  basis  elements,  instead  of 
sines  and  cosines  [8:2].  The  wavelets  must  be  oscillatory  (the  wave  part)  and  must  also 
decay  rapidly  (the  little  part).  An  example  of  such  a  wavelet,  called  the  Morlet  wavelet,  is 
shown  in  Figure  C.  1 . 


1 

'  1 

1 

r  ' 

Figure  C.  1 ;  Morlet  Wavelet 
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In  Fourier  analysis  the  pieces  used  to  decompose  and  reconstruct  a  signal  are  sines 
and  cosines  of  different  frequencies.  In  wavelet  analysis  the  building  blocks  are  translated 
and  dilated  (shifted  and  scaled)  versions  of  a  single  wavelet,  called  the  “mother  wavelet” 
[46:2],  Denoting  the  mother  wavelet  by  h(t),  the  rest  of  the  building  blocks  are  given  by 

where  a  is  the  scale  and  b  is  the  shift  [7:909],  The  mother  wavelet  h  must  satisfy  the 
condition 

r  h{x)dx  =  0  (C.2) 

J— 00 

A 

which  implies  that  /i(0)  =  0,  where  the  ^  denotes  Fourier  transform.  If  the  shift  h  and 
scale  a  vary  continuously,  the  continuous  wavelet  transform  based  on  this  mother  wavelet, 
h,  maps  a  one  dimensional  function  f(x)  to  a  two  dimensional  surface  Whf(a,b).  The 
mathematical  form  of  this  transform  is 

WJ{a,b)=r  f{x)h^^{x)dx  (C.3) 

J—oo  * 


If  the  shift  and  scale  are  restricted  to  a  discrete  sublattice,  the  continuous  wavelet 
transform  above  becomes  a  discrete  wavelet  transform.  The  wavelets  become  [7:910] 

hm,nit)  =  ^o”''^h{aQ’”t-nbQ)  m,n  gZ  (C.4) 

where 

a  =  ao 

b^nb^a^  (C.5) 

Z  is  the  set  of  integers 
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For  this  work  it  is  convenient  to  choose  ad=2,  and  hd=\,  although  other  choices  are 


possible.  With  this  choice  of  ao  and  bo  it  is  possible  to  construct  an  h(t)  such  that  hm.n 
form  an  orthonormal  basis  for  Z^(iR)  [7:911],  where  is  the  set  of  all  measurable 
functions  of  finite  energy,  defined  by 


^(9l)  =  |/ :/  is  measurable  and  |  |/(x)pfiSc  <  qo|  (C.6) 


Multiresolution  Analysis 

A  multiresolution  analysis  is  a  technique  of  representing  a  function  /eL^(91)  as  a 
limit  of  successive  approximations,  each  of  which  is  a  smoothed  version  of f  [7:913,  25], 
A  multiresolution  analysis  consists  of:  [7:915,  45:152-158] 

(i.)  A  group  of  embedded  subspaces  F^c  Z>^(5R),  such  as 

...cFj  eF,  cFq  cF_,  cF_2...  (C.7) 

(ii.)  These  subspaces  have  only  the  zero  element  in  common,  and  the  sum  of  the 
spaces  span  all  ofZ,^(9?): 


(C8) 

meZ  meZ 


(iii.)  If  a  function  is  an  element  of  one  space,  then  the  same  function  scaled  by  a 
factor  of  2  is  an  element  of  the  next  higher  space: 
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fix)eV„ 


fi2x)eV^_, 


(C.9) 


(iv.)  There  exists  a  function  Vo  such  that  shifted  versions  of  this  function  form  a 
basis  for  the  space  Vo: 


Vo  =span(^o,^,neZ} 


(CIO) 


The  first  property  above  means  that  any  feature  that  can  be  seen  at  coarse 
resolution  (like  Vj)  can  also  be  seen  at  finer  resolutions  (like  Vo).  Naturally,  the  finer 
resolutions  also  contain  more  detailed  features  as  well.  The  second  property  describes  the 
limits  of  this  ladder  of  subspaces.  At  the  coarsest  resolution  there  is  nothing  left  but  the 
zero  element  of  the  space,  while  in  the  limit  of  the  finest  resolution  the  Vm  converge  to 
Z^(9?).  The  third  property  is  the  same  scaling  feature  already  identified  as  a  property  of 
wavelets.  The  last  property  is  the  building  block  for  the  analysis.  For  a  fixed  dilation 
(scale),  translations  (shifts)  of  the  function  ^  are  basis  elements  for  that  subspace.  These 
functions  ^  are  called  scaling  functions,  or  father  wavelets.  Only  real  scaling  functions  will 
be  considered  here,  although  complex  ^  are  possible. 

Since  a  space  Vm  has  a  basis  given  by  then  any  element  of  the  space  can  be 
written  as  a  linear  combination  of  these  basis  elements.  That  is  if  x(t)  e  Vm  ,  then  there 
exist  coefficients  Cm,n  such  that 


neZ 


(C.ll) 
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The  Cm are  called  approximation  coefficients.  They  are  found  by  taking  the  inner  product 
of  the  function  x(t)  with  the  scaling  function 


c 


m^rt 


r  {t)dt 

J— oo  * 


(C.12) 


The  Dirac  bracket  notation  is  used  to  denote  the  inner  product  of  two  functions.  For 
,  a  projection  operator  Pm  can  be  defined  by 

=  (C.13) 

neZ 


which  gives  the  expression  for  that  part  of  f(t)  that  is  an  element  of  Since  all  the 
elements  in  ¥„  are  in  Vm-i  (by  the  construction  of  the  ladder  of  subspaces),  what  elements 
are  in  Vm-i  that  are  not  in  Fm?  Construct  a  projection  operator  Qm  that  identifies  such 
elements  by 


QJ  =  Pn.-xf-PJ  (C.14) 

By  the  Classical  Projection  Theorem,  the  0m/ is  orthogonal  to  Vm,  and  so  defines  a 
space  of  all  such  elements,  called  Wm,  whose  inner  product  with  elements  of  Vm  is  zero. 
In  set  notation,  that  is  [45:162] 


Wm={f  \  (/,v)  =  0VveFm}  (C.15) 


So  by  definition  of  Wm  and  0m 
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(C.16) 


Pm-J  =  PJ  +  Qmf 

C 


where  the  symbol  “©”  is  the  orthogonal  direct  sum  of  the  two  spaces.  The  direct  sum  is 
the  set  of  vectors  formed  by  adding  one  element  from  each  of  the  orthogonal  subspaces. 
Equation  C.16  is  a  key  point  in  the  multiresolution  analysis  theory.  Each  subspace  Vm-j 
consists  of  the  direct  sum  of  two  other  subspace,  Vm  and  This  idea  is  illustrated 
graphically  in  Figure  C.2. 


Figure  C.2:  Structure  of  Nested  Subspaces 

By  the  nature  of  the  ladder  of  subspaces  and  the  fact  that  W„  is  orthogonal  to 
the  Wm’s  are  all  mutually  orthogonal.  Furthermore,  the  combination  of  all  the  fVm’s  spans 
the  entire  space  of  functions: 


(C.17) 

/weZ 
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Now  identify  the  basis  for  each  JVm  as  {  neZ).  This  set  of  functions  forms  a  basis 
for  the  entire  space  Z,^(5R),  and  they  are  the  wavelets  we  seek.  The  details  of  the 
techniques  used  to  construct  the  scaling  function  and  the  wavelets  are  not 
necessary  to  understand  the  use  of  this  theory.  Several  different  types  of  scaling  functions 
and  wavelets  have  been  developed  in  the  literature.  These  include  the  Haar  basis  above, 
which  is  the  simplest,  as  well  as  Daubechies'  compactly  supported  wavelets  [7:980]. 
Illustrations  of  some  of  Daubechies*  wavelets  and  their  associated  scaling  functions  are 
shown  in  Figure  C.3.  The  different  wavelet-scaling  fimction  pairs  shown  are  labeled  by 
the  number  of  filter  coefficients,  or  “taps”,  needed  to  completely  specify  each  wavelet.  The 
functional  form  of  the  scaling  functions  and  wavelets  are  not  required  to  implement  the 
method  of  multiresolution  analysis,  as  will  be  seen  in  the  next  section. 
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Figure  C.3:  Illustration  of  Daubechies'  compactly  supported  wavelets  and  scaling  functions;  (a) 

2  tap  (Haar),  (b)  4  tap,  (c)  8  tap,  (d)  16  tap. 


Function  Decomposition  and  Reconstruction 

To  make  use  of  the  multiresolution  analysis  theory,  a  method  is  needed  to 
represent  the  function  to  be  analyzed  in  each  of  the  subspaces,  preferably  without 
performing  numerous  integrations.  To  develop  such  a  method,  define  the  coefficients  of 
expansion  for  each  level  “/«”  by 


(C.18) 


The  coefficients  Cm,„  are  found  by  taking  the  inner  product  of  the  projection  of f(t)  onto  Vm 
with  the  scaling  function  ^nt.n-  Applying  the  relationship  between  the  projection  operators 
in  Equation  C.  17  to  this  expression  for  c„,„  yields 


c  =(p  f ,6  \ 

m,n  \  mJ  >  r  m,n  f 


(C.19) 


However,  the  second  inner  product  is  zero  because  Qmf(t)  1  ^m.n  by  construction. 
This  leaves  a  recursive  relationship  between  the  coefficients; 

^m,n  ~  (C.20) 

<eZ 


Similarly  for  dm,n 
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(C.21) 


d„>,n={Qmf  ,¥m,n) 

,¥n,,n)-{Pmf  ^¥m,n) 

~  m-\,X  >  ¥  m,n  )• 

AeZ 

The  recursive  relations  above  show  all  that  is  needed  to  decompose  a  function  is  to 
find  the  inner  products  between  a  scaling  function  and  the  scaling  fimction  or  wavelet  of 
the  next  level  down.  Consider  the  explicit  form  of  one  of  these  inner  products 

- 1)  4(2- 1  -  n)dt  (C,22) 

Let  —  =  2“'"t-H,  Then 
2 

=  -2«])*  (C.23) 

Define  a  sequence  h(n)  by 


then 


h{n)  =  <l>{t-n)  dt 

(C.24) 

(C.25) 

Notice  this  relation  is  independent  of  the  level  ml  This  means  that  the 
coefficients  for  any  level  can  be  found  fi'om  the  coefficients  for  the  level  above.  (Also 
note  the  sequence  h(n)  defined  in  Equation  C.25  above  is  not  related  to  the  wavelet  h(t) 
defined  in  the  previous  section.  This  conflicting  notation  is  standard  within  the  wavelet 
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literature.)  The  coefficients  required  to  specify  the  Daubechies'  wavelets  shown  in  Figure 
C.3  are  listed  in  Table  C-1. 


Table  C-1 :  Daubechies'  wavelet  coefficients  [7:980] 


h-2:  0.7071067811865 

0.7071067811865 


h=4:  0.482962913145 

0.836516303738 
0.224143868042 
-0.129409522551 


h=8:  0.230377813309 

0.714846570553 
0.630880767930 
-0.027983769417 
-0.187034811719 
0.030841381836 
0.032883011667 
-0.010597401785 


h=16:  0.054415842243 

0.312871590914 
0.675630736297 
0.585354683654 
-0.015829105256 
-0.284015542962 
0.000472484574 
0.128747426620 
-0.017369301002 
-0.044088253931 
0.013981027917 
0.008746094047 
-0.004870352993 
-0.000391740373 
0.000675449406 
-0.000117476784 
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Thus,  there  is  a  recursion  relation  between  levels  that  can  be  exploited  to  calculate 
the  coefficients: 

o^„='Lc^A(-2n)  (C.26) 

<gZ 

A  similar  relation  exists  for  the  detail  coefficients,  dm_„\ 

g{n)  =  2“’^' J  ^  ^(i)  (f>{t  -  n)  dt  (C,27) 

The  decomposition  algorithm  in  this  multiresolution  analysis  requires  only  the 
filters  h  and  g  to  define  the  wavelets  to  be  used.  In  fact,  the  g  filter  can  be  derived  from 
the  h  filter,  so  only  one  sequence  is  needed  to  completely  determine  the  wavelets.  This 
relationship  between  h  and  g  is 

g{n)  =  {-\fh{n-\)  (C.28) 

Using  the  two  recursion  relations  for  c„,,„  and  any  function  can  be 
decomposed  into  its  approximation  and  detail  coefficients  from  an  initial  sequence  Co. 
Figure  C.4  illustrates  the  relationship  between  these  coefficients.  Note  the  structure  of  the 
coefficients  is  the  same  as  the  structure  of  the  underlying  spaces,  Vm  and  Wm. 

Now  that  the  technique  used  to  decompose  a  function  has  been  identified,  how  is 
the  original  function  reconstructed  from  these  coefficients?  As  can  be  seen  from  Figure 
C.4,  only  the  coarsest  (lowest)  approximation  coefficient  and  the  detail  coefficients  are 
needed  to  reconstruct  the  original  Co. 
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The  recursive  relation  for  reconstruction  is 

=  S  -  2")  +  Z  -  2«)  (C.29) 


A  further  simplification  to  the  decomposition-reconstruction  algorithm  can  be 
made  by  treating  the  initial  sequence  of  coefficients  as  part  of  a  periodic  sequence. 
Wavelet  transforms  have  the  property  that  the  transform  of  a  periodic  function  is  also 
periodic.  In  the  case  of  the  multiresolution  decomposition,  if  the  top  sequence  co  has  N 
values,  and  A  is  a  power  of  2  (i.e.,  A=2^  for  some  integer  M),  then  the  next  level  coarse 
coefficients  c;  will  be  2^  *  periodic,  as  will  the  details  for  that  level,  dj.  So  the  lowest  level 
of  decomposition  consists  of  a  1 -periodic  coarse  coefficient  sequence  (a  constant),  and  a 
1 -periodic  detail  coefficient  sequence. 
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